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Variational Perturbation Theory for Fokker-Planck Equation with Nonlinear Drift 
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We develop a recursive method for perturbative solutions of the Fokker-Planck equation with nonlinear drift. 
The series expansion of the time-dependent probability density in terms of powers of the coupling constant is 
obtained by solving a set of first-order linear ordinary differential equations. Resumming the series in the spirit 
of variational perturbation theory we are able to determine the probability density for all values of the coupling 
constant. Comparison with numerical results shows exponential convergence with increasing order. 
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I. INTRODUCTION 

In many problems of physics, chemistry and biology one 
has to deal with vast numbers of influences which are not 
fully known and are thus modelled by noise or fluctuations 
IHSbfiU- The stochastic approach to such systems identi- 
fies some relevant macroscopic property x governed by a drift 
coefficient K{x), while the microscopic degrees are averaged 
to noise, entering as a stochastic force. In the case of addi- 
tive noise this simply leads to a diffusion constant D. Thus 
we arrive at the Fokker-Planck equation (FP) |3l for the den- 
sity IP (x, f ) of the probability to find the system in a state with 
property x at time t: 

j^'Pi^.t) = -|- [K{x)T{x,t)] +D^¥{x,t) . (!) 

Once £P {x,t) is known from solving we can calculate the 
ensemble average of any function O (x) according to 



(0(x(f))>=/ 0{x{t))¥{x,t)dx. 



(2) 



In this paper we consider a stochastic model with the nonlin- 
ear drift coefficient 



K{x) — —yx — gx^ 



(3) 



Such models are studied, for instance, in semiclassical 
treatments of a laser near to its instability threshold |5, 6], 
where the variable x is taken to be the electric field. The 
damping constant y is set proportional to the difference 
between the pump parameter o and its threshold value Othr, so 
that the laser instability corresponds to y = 0. The coupling 
constant g > describes the interaction between light and 
matter within the dipole approximation, and the diffusion 
constant D in the FP equation characterizes the spontaneous 
emission of radiation. While the linear case with g = 0, 
corresponding to Brownian motion with damping constant y. 
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can be solved analytically, there is no exact solution for the 
nonlinear case with g > 0. But we can find a solution in form 
of a series expansion of the probability density iP(x,f) in 
powers of the coupling constant g. This series is asymptotic, 
i.e. its expansion coefficients increase factorially with the 
perturbative oder and alternate in sign. 

Such divergent weak-coupling series are known from 
various fields of physics, e.g. quantum statistics or crit- 
ical phenomena, and resummation techniques have been 
invented to extract meaningful information in such situations. 
Powerful tools among them are the variational methods, in- 
dependently studied by many groups, (see e.g. the references 
in tZil). A simple example is the so-called 5-expansion of 
the anharmonic quantum oscillator, where the trial frequency 
of an artificial harmonic oscillator introduced to maximally 
counterbalance the nonlinear term, is optimized following 
the principle of minimal sensitivity 0]. It turns out that the 
5-expansion procedure corresponds to a systematic extension 
of a variational approach in quantum statisti cs P. fiol fill fl^ 
to arbitrary orders as developed by Kleinert \V, 

mm 

, now 

being called variational perturbation theory (VPT). In recent 
years, VPT has been extended in a simple but essential way 
to also allow for the resummation of divergent perturbation 
expansions which arise from renormalizing the (|)'*-theory of 
critical phenomena QUI [13 El. The most important 
new feature of this field-theoretic variational perturbation 
theory is that it accounts for the anomalous power approach 
to the strong-coupling limit which the 5-expansion cannot do. 
This approach is governed by an irrational critical exponent 
as was first shown by Wegner |19] in the context of critical 
phenomena. In contrast to the 5-expansion, the field-theoretic 
variational perturbation expansions cannot be derived from 
adding and subtracting a harmonic term. Instead, a self- 
consistent procedure is set up to determine this irrational 
critical Wegner exponent. The theoretical results of the 
field-theoretic variational perturbation theory are in excellent 
agreement with the only experimental value available so far 
with appropriate accuracy, the critical exponent a governing 
the behaviour of the specific heat near the superfluid phase 
transition of ^He which was measured in a satellite orbiting 
around the earth 

Recently, the VPT techniques have been applied to Markov 
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theory, approximating a nonlinear stochastic process by an 
effective Brownian motion I23I1 . This is achieved by 
adding and subtracting a hnear term to the nonhnear drift 
coefficient (|3}, where the new damping constant is taken as 
the variational parameter. In the present paper we extend 
this result to higher orders, which we have made accessible 
by our recursive approach l24ll . By doing so, we are able to 
show that VPT makes Markov theory converge exponentially 
with respect to order, a phenomenon known for various other 
systems EUlEllllllllSIIlll. 

The paper is structured as follows. In Sec. HI] we review 
some properties of the FP equation which are essential for our 
discussion. In Sec.|n]we present the asymptotic perturbation 
expansion for the normalization constant of the stationary so- 
lution of the FP equation with the drift coefficient (|3ji and its 
variational resummation. This simple case already illustrates 
in an introductory fashion all features of the upcoming treat- 
ment of the time-dependent problem. In Sec. llVl we perturba- 
tively solve the FP equation with a nonlinear drift coefficient 
(|3} for Gaussian initial distributions by means of a double ex- 
pansion with respect to the coupling strength g and the vari- 
able X. In Sec. |V] we apply VPT to the resulting divergent 
weak-coupling expansion of the probability density to render 
the results convergent for all values of the coupling constant g. 
Furthermore, we discuss the exponential convergence of our 
variational method with respect to the order The paper closes 
with a summary in Sec lVll 



with probability current 

Six,t) =K{x)¥{x,t)-D^¥{x,t). (7) 
ax 

For natural boundary conditions, where S{x,t) ^ as x ^ 
±00, probability conservation is thus guaranteed: 



— J dx¥{x,t)=0. 



(8) 



In the case of additive noise, the drift coefficient K{x) can be 
derived from a potential 



<P{x) 



K{y)dy 



such that 



K{x) ^ -D^'{x). 



(9) 



(10) 



If the potential satisfies 't'(x) +°o for x — > ±00, the proba- 
bility density £P {x,t) approaches a stationary values iPstat(-«) in 
the long-time limit, i.e. 



i'statW = lim2'(x,f) =Ne 



where the normalization constant N is found to be 

-1 

i / " / N I 

exp 



^ / K{y)dy 



dx 



(11) 



(12) 



II. FOKKER-PLANCK EQUATION 



B. Brownian Motion 



In this section we fix our notation by reviewing the main 
properties of one-dimensional stochastic Markov processes. 



A. Definitions 

A stochastic Markov process x{t) is described by an ordi- 
nary stochastic differential equation which is of the Langevin 
form 



x{t) ^ a{x{t)) +b{x{t))T{t) , 



(4) 



where r(f ) is a Gaussian distributed stochastic force with zero 
mean and 5-correlation: 



(r(f))=0, (r(r)r(r'))=28(r-f' 



(5) 



For the special case of the coefficient b{x) = b being constant, 
the stochastic system is said to be driven by additive noise, 
and the time evolution of its probability density 2'(x,f) is de- 
scribed by the Fokker-Planck equation Q with drift coeffi- 
cient K{x) = a{x) and diffusion coefficient D = b^ 0]. The 
FP equation Q has the form of a continuity equation 



— 'P{x,t) + —S{x,t) =0 
at ax 



(6) 



For Brownian motion, defined by a linear drift coefficient 
K{x) = -yx (13) 
with damping y > 0, we have a harmonic potential: 

(14) 



4>(x) = ^x^ 
^ ' ID 



and the FP equation Q has a solution in closed form. With 
initial condition 



£P (x,xo;0) = 5(x — xo) 



(15) 



the solution is 



iP(x,x';f) = ,/ — 



X exp 



(x-x!e-^'f • 
"2D(l-e-2Yf)/y 



(16) 



It represents the Green function of the FP equation with the 
drift coefficient (I13> . describing the evolution of the proba- 
bility density starting with an arbitrary initial distribution ac- 
cording to 



£P(x,f)= y 'P{x,x';t-t')¥{x' ,t')dx' . 



(17) 
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For instance, the initial Gaussian probability density 



has the time evolution 



1 



exp 



2o2 



(18) 



27t[D-(D-yo2)e-2T] 



X exp 



(19) 



2D-(D-Yo2)e-27'_ 
In the long-time limit it approaches the stationary distribution 

fPstat(x) = Um ¥y{x,t) = xf^^W f-Sr) (20) 



which turns out to be independent of the parameters fj and O 
of the initial distribution ( I18> . 



C. Nonlinear Drift 



Consider now an additional nonlinear drift Q so that the 
potential (|9} becomes 



4>(x) 



— X H X 

2D 4D 



(21) 



While there exists no closed solution of the corresponding FP 
equation Q 

|:2'(jc,0 = ^ [{yx+gx^)¥{x,t)] +D^(p{x,t) 

its stationary distribution il 1> is given by 

2'statW=N(^)expf--^x2-A 



(22) 



2D 



AD 



with normalization constant 



(23) 



(24) 



Performing the integral in (I24> . we obtain 

^exp(-5^) /Ki/4(5gj 
uAl/4 



2^ 



1/4 (8Dg )+'l/4U£lg ) 



;y>0 
; Y=0 

;y<0, 



(25) 



where Iv(z) and Kv(z) denote modified Bessel functions of 
V-th order of first and second kind, respectively ll32ll . 



III. NORMALIZATION CONSTANT 

The diverging behaviour of a perturbation series as well as 
the method of VPT to overcome this problem can already be 
studied by considering the normalization constant \2Al of the 
stationary solution \23\ . To simplify our discussion we as- 
sume in this section without loss of generality D =\. 



n 


V^27tfl„ 







1 


0.390062251089407 =r(3/4)/7t 


1 


3/4 


0.131836797004050253244 


2 


-87/32 


-0.004 1 98 378 378 722 963 623 


3 


2889/128 


-0.001419006213792844574 


4 


-581157/2048 


0.000536178450689882683 


5 


38668509/8192 


-0.000093 437 5 1 1 028 762 876 



TABLE 1: The first 6 coefficients of the weak- and strong-coupling 
expansion <26> and <29> . respectively, for Z) = y = 1. 



A. Weak- and Strong-Coupling Expansion 

The weak-coupling expansion 



N, 



N 



weakfe) = L 

«=0 



(26) 



follows from ( I24> by expanding the exp(— ^x^/4)-term in the 
integrand; 



/gf (-g)^r(l/2 + 2^) 
to k\ fl^+'-'^ 



-1 



(27) 



On the other hand, a rescaling x {'^/gY^'^y in the normal- 
ization integral J24> leads to 



1/4 



(28) 



so the strong-coupling expansion 



N, 



(M) 
Strong 



M 



,-'n/2 



(29) 



follows from ( I28> by expanding the term exp(— yy2/g'/2) in 
the integrand 



^^Ug) = (I 



1/4 



^ (-#£(1/4 + ^2) 



H -1 



2k\ 



(30) 



Note that the weak-coupling expansion ( I26> contains integer 
powers of g, whereas dimensional arguments lead to rational 
powers of g for the strong-coupling expansion \29\ . The first 
six coefficients an,b,„ are given in Table U and the respective 
expansions ( I26l l and ( I29> are depicted in Figure^ While the 
weak-coupling expansion provides good results for small val- 
ues of g, the strong-coupling expansion describes the behavior 
for large values of g. For intermediate values of g both series 
yield poor results. 
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FIG. 1: Weak- and strong-coupling expansions 1261 and <29> up to 
tlie 5th order, respectively, as well as the exact normalization constant 
l24l(«««)for£) = Y= 1. 



B. Variational Perturbation Theory 

Despite its diverging nature, all information on the analytic 
function N{g) is already contained in the weak-coupling ex- 
pansion i26\ . One way to extract this information and use 
it to render the series convergent for any value of the cou- 
pliiig co nstant g is provided by VPT as developed by Kleinert 
lll4l [Tsl llSll . This method is based on introducing a dummy 
variational parameter on which the full perturbation expansion 
does not depend, while the truncated expansion does. The op- 
timal variational parameter is then selected by invoking the 
principle of minimal sensitivity 0], requiring the quantity of 
interest to be stationary with respect to the variational param- 
eter In our context l22ll23L 12411 . this dummy variational pa- 
rameter can be thought of as the damping constant K of a trial 
Brownian motion with a harmonic potential Kx^/2, which is 
tuned in such a way, that it effectively compensates the nonlin- 
ear potential. In order to introduce the variational parameter 
K, we add the harmonic potential of the trial Brownian 

motion to the nonlinear potential i2H and subtract it again: 



(31) 



By doing so, we consider the harmonic potential Kx^/2 as the 
unperturbed term and treat all remaining potential terms in 
OH as a perturbation. Such a formal perturbation expansion 
is performed for the normalization constant (I24> according to 



N(^) 




exp 



-5 



x" 



dx 



8=1 



where the additional parameter 8 is introduced in the spirit of 
the 5-expansion (see, for instance, the references in 0]). Due 
to the rescaling x ^/P(8) with the scaling factor 



P(5) = V[K + 5(y-K)]/Y, 



(33) 



the weak-coupling expansion (I26> of \2A\ leads to: 



N(g)=P(8)^a„ 

n=0 



(34) 



8=1 



Expanding and truncating the series at order in 8, we obtain 
the Mh variational approximation to the normalization con- 
stant N(g): 



\/2-2n 




(35) 



Equivalently, the variational expression ( I35> also follows di- 
rectly from the weak-coupling expansion \2bl . To this end 
we remark that, due to dimensional arguments, the respective 
coefficients depend via 



(36) 



on y where a,, denote dimensionless quantities. Treating in 
(13 H the harmonic potential Kx^ /2 as the unperturbed term and 
all remaining potential terms as a perturbation, corresponds 
then to the substitution 



Y^K(l+§r) 



with the abbreviation 



(37) 



(38) 



Thus inserting J37> in the weak-coupling expansion J26t . 
( I36l l. a reexpansion in powers of g up to order together with 
the resubstitution ( I38> leads to a rederivation of (I35> . 

As the normalization constant N(^) does not depend on the 
variational parameter K, it is reasonable to ask for the trun- 
cated series J35> to depend on K as little as possible. There- 
fore, we invoke the principle of minimal sensitivity 0] and 
demand 



0. 



(39) 



The optimal variational parameter K^o2{g) leads to the Mh 



opt 

order variational approximation ^''y^{gMo^i{s)) to the nor- 
malization constant N(g). In case that ( I39> is not solvable, we 

determine Kq^j is) from the zero of the second derivative in 
(32) accordance with the principle of minimal sensitivity |3l: 



33yNfp'T(g,K) 



0. 



(40) 



If it happens that ( I39l l or \AQ\ have more than one solution, 
we select that particular one which is closest to the solution 



5 




0.5 



1.5 



FIG. 2: First-oder variational result ( — ) compared with exact nor- 
malization constant (ooo). First-order weak- and strong-coupling 
approximations are shown by dashed ( — ) and dotted (■ ■ ■) lines, re- 
spectively. 



in the previous variational order. 

As an illustrative example we treat explicitly the first varia- 
tional order where we obtain 



(1) _ 2K(K + y) + 3g 

so that the first derivative with respect to K 

has the two zeros 



(41) 



K 



opt ig)^^{j±\/^8g + f 



(43) 



Since the variational parameter has to approach y for a van- 
ishing coupling constant g, we select from ( I43> the solution 

K^ipf 'te)- The resulting optimized result N^'^T(^'i^ipt^'(^)) 
is shown in Figure |2] We observe that in this parame- 
ter range already the first-order variational approximation 



NypT fe' ^apt^' (§) ) is indistinguishable from the exact normal- 
ization constant N{g). 



C. Exponential Convergence 

In order to quantify the accuracy of the variational approx- 
imations, we study now, in particular, the strong-coupling 
regime g ^ In first order, the insertion of ( I43> in ( I41> 
leads to the strong-coupling expansion i29\ with the leading 
coefficient 



yields the leading strong-coupling coefficient within an 
accuracy of less than 1 %. 

To obtain higher-order variational results for this strong- 
coupling coefficient bo, we proceed as follows. From the first- 
order approximation ( I43> we see that the variational parameter 
has a strong-coupling expansion K^pf^' (g) = y/g3V2/2 + 
whose form turns out to be valid also for the orders N > I. 



Inserting the ansatz K^pj (g) — Kq \/g + ... in ( l35t , we obtain 
the Mh order approximation for the leading strong-coupling 
coefficient bo'. 



{N) \ 1/2- 2« 




(45) 



The inner sum can be performed explicitly by using Eq. 



(0.151.4) in Ref. |3 



4'"(Kr) = L(-ir" 

n=0 



-1/2 -2n 
N-n 



■ (46) 



^0 

y 



In order to optimize (146 > we look again for an extremum 



= 



(42) or for a saddle point 



(47) 



(48) 



It turns out that extrema exist for odd orders A^, whereas even 
orders lead to saddle points. The points of Figure |3] show 

the logarithmic plot of the relative error |^o''''(Kq— — bo\/bo 
when using the smallest zero of the first and second deriva- 
tive, i.e. iA-ll and ( I48t . respectively. We observe that the rel- 
ative error depends linearily on N^^^ up to the order = 120 
according to 



bo 



(49) 



where the fit to the straight line —a — PA^'^^ leads to the 
quantities a = 0.139714 and p = 1.33218. Thus we have 
demonstrated that the variational approximations for the lead- 
ing strong-coupling coefficient converge exponentially fast. 
Note that the speed of convergence is faster than the exponen- 
tial convergence of the variational results for the ground-state 
energy of the anharmonic oscillator l2^E3l . 



(1) 



2 

97t2 



1/4 



: 0.387. 



(44) 



Comparing ^ with bo = r(3/4)/7t w 0.390 (see Table Q, 
we conclude that first-order variational perturbation theory 



IV. RECURSION RELATIONS 

Now we elaborate the perturbative solution of the FP equa- 
tion i22\ with the initial distribution (I18t . By doing so, we 
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-0= 




10 11 



FIG. 3: The points show the logarithmic plot of the relative error 
I ^Q^' (Kq"' ) ~ ^0 1 /bo when using the smallest zero of the first (•) and 
the second derivative (o), respectively, against W'/^. The solid line 
represents a fit to the straight line —a — PA''/^. 



where we have introduced the abbreviation 



(54) 



For a coupling constant § > 0, we solve ( 15 U by the ansatz 

ip(x,t) = 2'y(x,t)^(x,t), (55) 

so that the remainder q{x,%) fulfills the partial differential 
equation 



3gx 



2 , Jgx^ -Jg'^XQX^ 



D(t2-1) 



Y(t2 + i)x-2yt:xo , ^ 
+^-" 



—q{x,z)+D^q{x,%). (56) 



Then we solve ( I56> by expanding ^(x,t) in a Taylor series 
with respect to the coupling constant g, i.e. 



follow the notion of Ref. 113011 and generalize the recursive 
Bender- Wu solution method for the Schrodinger equation of 
the anharmonic oscillator l33ll . thus obtaining a recursive set 
of first-order ordinary differential equations 12411 . 



A. Time Transformation 

At first we perform a suitable time transformation which 
simplifies the following calculations: 



q{x^) = Y^g"qn{x^), 

n=Q 



(57) 



where we set qo{x,z) = 1. Thus the expansion coefficients 
qnixjX) obey the partial differential equations 



-y^^qnix^z) 



3x- 



2 yx^ — yixQX^ d 



D(t2-1) 



^„_i(x,t) 



(50) 



Thus the new time x runs from Tq to when the physical time 
t evolves from to °o. Due to i5Q\ the FP equation J22> is 
transformed to 

-Yi|^£P(x,t) = |- [{yx + gx^)'P{x,z)] +D^¥{x,t:). (51) 
Furthermore, the initial distribution ( I18t reads then 



rp(x,T = To) = 



1 



:exp 



2a2 



(52) 



and contains P{x,z = 1) = 5(x — /j) as a special case in the 
limit a2 0. 



B. Expansion in Powers of g 

If the coupling constant g vanishes, the solution of the ini- 
tial value problem i5ll . ( I52t follows from applying the time 
transformation ( I50> to ( I19> , i.e. 



2'y(a:,t) 



27tD(l-T2) 



exp 



Y (x-xot)2 



2D 1 - t2 



(53) 



C. Expansion in Powers of jc 

It turns out that the partial differential equations (I58> are 
solved by expansion coefficients qn {x, x) which are finite poly- 
nomials in x: 



M„ 



qn{x,l) ^ a„,i(T)x^ 



(59) 



k=0 



Indeed, inserting the decomposition (I59> in \5%\ . we deduce 
that the highest polynomial degree is given by M„ = An. Note 
that in case of /j = the partial differential equations ( I58> are 
symmetric with respect to x — x, so that only even powers x 
appear in ( I59> . From ( I58> and ( I59> follows that the functions 
(X„,k{z) are determined from 

a , , k(z^ + l) 

— a„.A.(T + — ^ 7T0C„,i T = r„.i(T) , (60) 

OT T(T^ — 1 j 

where the inhomogeneity is given by 

k+l D{k + 2)ik+l) 

'"n.A'('C) = —C^n-Lk-li"^) — a«.i+2('C) 



J^0a«-1,/S:-3('C) 0C„_i,t._4(T) 2xo(fc+l) 



DrT2_l) 



D(t3-t) 



t2-1 



an,/i+i(i:)-(61) 
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k=0 
n=0 -4- 



n=2 @^ 



k=l 



k=2 k=3 k=4 k=5 



k=6 k=7 



k=8 




FIG. 4: Recursive calculation of the functions a„ i^{z): ao o{l) = 1 
is the only coefficient which is a priori nonzero (0). For each n 
the coefficients are successively determined for k = 4n, . . . ,0 (|0|). 
Each step necessitates (■ " " ») only those coefficients which are already 
known (®). 



D. Recursive Solution 

The Eqs. ( I60t and i6H represent a system of first-order 
ordinary differential equations which can be recursively 
solved according to Figure |3 By doing so, one has to take 
into account OC„.a:(t) = if « < or ^ < or ^ > 4«. The 
respective positions in the (n,fc)-grid of Figurel^are empty or 
lie outside. Iteratively calculating the coefficients a„^k{z) in 
the nth order, we have to start with k = An and decrease k up 
to A: = 0. In each iterative step the inhomogeneity (I61> of ( I60> 
contains only those coefficients which are already known. 
The only coefficient which is a priori nonzero is ao.o (x) = 1 . 

Applying the method of varying constants, the inhomoge- 
neous differential equation j60> is solved by the ansatz 



that A„ ,t(To = 1) = 0, otherwise (Xn^ix) in \62\ would posses 
for T = 1 a pole of order k. For A: = this argument is not 
valid as the fraction in \62\ is then no longer present. The in- 
tegration constant A„. 0(^0 = 1) follows from considering the 
normalization integral of \55\ for x = 1 together with \51l . and 
(El, i.e. 



4n 



(66) 



n=\ k=0 

Indeed, taking into account (|62} for k = Q leads to 

An.oiT: = 1) = a„.Q{^ = 1) = - E ^"-ki^ = l)''^o- (67) 

k=l 

Thus for k y^O the coefficients a„.A (T) are still given by il}, 
whereas for A: = we obtain 



rz 4)1 

a„,o(T:) = / r„.o{a)da~ Y a„.k{i = 1)4- (68) 
•^1 k=\ 



E. Cumulant Expansion 

The weak-coupling expansion i55\ . i51\ of the probability 
density iP(jc,T) has the disadvantage that its truncation to a 
certain order could lead to negative values. To avoid this, 
we rewrite the weak-coupling expansion ( I55t . ( I57t in form of 
the cumulant 



!P(x,t) = e 



p(A-,0 



(69) 



where the exponent p{x,z) is expanded in powers of the cou- 
pling constant g: 



where A„ ^(t) turns out to be 



(t2- 



(62) 



(63) 



The integration constant A„ j.(to) is fixed by considering the 
initial distribution From O, (EJ, (EJ, and (|59j fol- 
lows 



4n 



i = i + E^"E«"i(^o)/, 

n=l k=Q 



(64) 



SO we conclude a„,i(To) = and thus A,2,i:(To) — for « > 1. 
Therefore, we obtain the final result 



an,ki^) = 



r„,k{a)o ''{a- -If da, 



(t2-1)^..„ 
where r„^ii{a) is given by ( l6ll 



(65) 



Note that the special case ^ with the initial distribu- 
tion P{x,x = 1) = 5(x — ju) has to be discussed separately, as 
then (I64> is not valid. In this case we still conclude for k^Q 



p{x, t) = In 2'y(jc, x) + £ g" Pn (J^, t) . (70) 

n=\ 

The respective coefficients p„{x,z) follow from reexpanding 
the weak-coupling expansion ( I55> . according to ( I69> , 
(I70> . However, it is also possible to derive a recursive set of 
ordinary differential equations whose solution directly leads 
to the cumulant expansion ( l69l , dTOl . To this end we pro- 
ceed in a similar way as in case of the derivation of the weak- 
coupUng expansion and perform the ansatz 



2n+22n+2-k 
Pn{x,^) = E E KkMx'^x'q. 

k=Q 1=0 



(71) 



The respective expansion coefficients ^„,k.i{i) follow from a 
similar formula than ( I65> , i.e. 

= , / /\„.A,,(o)o-*(o'- l)'£/a, (72) 

(Z - I) Jxo 

where the functions s„_k.i{'^) are given for n = 1 by 

, 3 8(t.35/.i 5a:.45/,o 

^1,.,/ W = --8a.,25.o + - 

, 2(A:+1) „ D(fc + 2)(fc+l) „ 

— rPi,-t+i,/-i(T:) — Pi,-t+2,/(i: , (73) 
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and for « > 2 by Iterating ( I72t -( l74t one has to take into account that ^n^k.i (t) 

vanishes if one of the following conditions is fulfilled: « < 0; 
, 2{k+l) k <Q OT k>2n + 2; I <0 OT I >2n + 2-k- k + l odd. By 

yi T-^ — 1 inverting the time transformation (15 Oil, the expansion coeffi- 



yi ' ' 'V- — \ ' ' inverting the time transformation ( 15 Ot . the expansion coeffi 

D{k + 2)(k+ 1) D cients are finally determined as functions of the physical time 

-Pn,jt+2./('c) H 1^ ^ y + 2) f • For the first order n~\ one finds the following expansion 



'^m=\j=i coefficients Pi, <../(f) 



X L ^n^.j.i{'^)^n-m.k-i+2.l-M ■ (74) 
i=0 



Pi,o,o(0 
PlA2(f) 
Pi, 0,4 (f) 
Pl,l,lW 
Pl,l,3(f) 
Pl,2,0(f) 
Pi, 2,2 (f) 
Pl,3,l(f) 
Pl,4,0(f) 



3Z)(g-(-5+f-"^-4>7+c-'T(4-8fY))+8DY(l+fY+e^'i'(-l+fY))c^-2f (l-e-'^+2fY)a'') 
4f (D(-l+i?2'i')+Yo2)- 

(0( 1 +e*'''(-5+Ati) +e-'T(4+8fY)) -2y( 1 -e'"^+4e-'i'fY)a- ) 
2Y(£>(-l+e2'T)+Ya2)' 

( 1 -6e^'T+2g'"^+f ■"^(3 - 1 2fY) ) 
4(D(-l+e-'T)+ya-)"' 

-3De'T'(20^(4+3fY+e'''^(-2+fY)+g-'^(-2+8fY))-OY(13-c''''y+12fY+4e^'^(-3+4rY))g-+3f (l-c^'^+2fY)g*) 

2Y(I>(-l+e-''')+ya2)-' 

-(D^e'T(a(-l+c''"'+e-"^(9-12fY)-3e^'^(3+4fY))+3Y(l-e'"^+4e-'^fY)g^)) 
2{D(-l+c2'T)+Yo2)'' 

3(D-^(l+e*'T'(-5+4fY)+g-'^(4+8>y))-D-Y(3+e-''T(-7+4>7)+4e^'T(l+4fY))g-+gf (3-c'''T+e-'T 

2Y(D(-l+e'2'T)+Yo2)-' 

-3De-'^(D-(3+2>y+8c-'^fy+f'''i'(-3+2fy))-DY(5-e*'V4>y+e-'i'(-4+8>y))a^+r(l-e-'^+2fY)a'') 

e'^(g-^(l-e'''^+3f-''T'(-3+4fY)+3e^''y(3+4fy))-30^y(l+e''''y(-5+4fY)+e^'y+8ty))a-+ 

l(D{-\+e^'1)+yc-)'' 

-(e'-'\D\2-6e*'-l+e^'-<+ie-'\l+An))-6D^-l(\-e*'-<+Ae'-'-<ti)o^+bDf^^^ 

4(D(-l+e2'T)+ya2)'' 

TABLE II: First order cumulant expansion coefficients Pi,/t,/(f). 



They are plotted for y > in Figure |5l where we distinguish 
three time regimes from their qualitative behavior. In the 
limit f all expansion coefficients ^\,k.i{t) vanish as al- 
ready the probability density of the Brownian motion M9\ in 
\69\- M\\ leads to the correct initial distribution ( I18> . In the 
opposite limit t ^ oa the only nonvanishing expansion coeffi- 
cients Pi, <:,/(?) read 

3D I 

Y>0: limpi,o,o(r) = -3, Um pi,4,o(f) = -— , (75) 

f^oo (^oo 4£) 



so that (I69>-(I71> reproduces the correct stationary solution 
( I23> . ( I26> up to first order in g. For intermediate times t we 
observe that all expansion coefficients pi, <:,/(?) show a non- 
trivial time dependence. Note that the number of coefficients 
Pn.*:,/(0 which have to be calculated in the «th order is given 
by L^io(2^+ 1) = (n +2)^, thus it increases quadratically. 
The expansion coefficients P„,<:,/(f) up the 7th order can be 
found in Ref. 13. 



I 

0.8 




1 2 3 4 J 5 

FIG. 5: Time evolution of the expansion coefficients Pi,i-,;(0 from 
Table II for £) = 1,Y= 1, = and o = 0.5. 
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V. VARIATIONAL PERTURBATION THEORY 



also to /7'^)(x,f;K) 



In this section we follow Refs. I22ll24ll and perform a vari- 
ational resummation of the cumulant expansion in close anal- 
ogy to Section ITllBI By doing so, we variationally calculate 
the probability density fP (x, f ) for an arbitrarily large coupling 
constant g with y > (anharmonic oscillator) and y < (dou- 
ble well). In both cases, we obtain probability densities, which 
originally peaked at the origin, turn into their respective sta- 
tionary solutions in the long-time limit. 



A. Resummation Procedure 

We aim at approximating the nonlinear drift coefficient (|3j 
by the linear one — Ki of a trial Brownian motion with a damp- 
ing coefficient K which we regard as our variational parameter 
To this end we add —Kx to the nonlinear drift coefficient (|3j 
and subtract it again: 



K{x) ~ — Kx — gx^ — (y — k) 



(76) 



By doing so, we consider the linear term —Kx as the unper- 
turbed system and treat all remaining terms in ( I76l l as a per- 
turbation. Such a formal perturbation expansion is performed 
by introducing an artificial parameter 5 which is later on fixed 
by the condition 8=1: 



K{x) = -KX - 8 [gx^ + (y- k)x] 
Performing the rescaling 



P'(8)l 



P(8)' 



p2(8) ' 



(77) 



(78) 



with the scaling factor ( I33t . the FP equation Q with the drift 
coefficient ilH is transformed to the original one i22\ . Due 
to dimensional reasons also the parameter /j,o of the initial 
distribution have to be rescaled according to 



(79) 



m 



m 



The rescaling ( I78> , ( I79> is applied to the cumulant expansion 
( I69> . (^Oj. After expanding in powers of 8 and truncating 
at order A^, we finally set 8 = 1 and obtain some Mh order 
approximant p'^' {x,t;K) for the cumulant. 

Equivalently, the same result follows also from the weak- 
coupling expansion of the cumulant ( I69l l, ( 17 Ot : 



(x,f) = ln£Py(x,f) + £ pnix,t)g" . (80) 
,1=1 



Treating in (I76> the linear drift coefficient —Kx as the 
unperturbed term and all remaining terms as a perturbation 
corresponds then to the substitution (|^} with the abbreviation 
( I38> . Thus inserting (|^} in ( I80> . a reexpansion in powers of 
g up to the order together with the resubstitution ( I38> leads 



If we could have performed this procedure up to infinite 
order, the variational parameter K would have dropped out of 
the expression, as the original stochastic model (|3} does not 
depend on K. However, as our calculation is limited to a finite 
order A^, we obtain an artificial dependence on K, i.e. some A^th 
order approximant p*'^'(x,r;K), which has to be minimized 
according to the principle of minimal sensitivity Q]. Thus 
we search for local extrema of p^^'>{x,t;K) with respect to K, 
i.e. from the condition 



8p(^)(x,f;K) 



dK 



(81) 



K=*:iB,'(A-.0 



It may happen that this equation is not solvable within a cer- 
tain region of the parameters x,t. In this case we look for zeros 
of higher derivatives instead in accordance with the principle 
of minimal sensitivity 0], i.e. we determine the variational 
parameter K from solving 



dK" 



(82) 



■^opt.mV 



The solution K^} {x,t) from (18 U or ( I82> yields the variational 
result 



P{x,t) 



exp 




x,t;Kip}{x,t)^ 






-piN) 


(^x ,f;KQp{'(jc ,f)^ 


)] 


dx' 



(83) 



for the probability density. Note that variational perturbation 
theory does not preserve the normaUzation of the probability 
density. Although the perturbative result is still normalized 
in the usual sense to the respective perturbative order in the 
coupling constant g, this normalization is spoilt by choosing 

an x-dependent damping constant Kg'^t (x,f). Thus we have 
to normalize the variational probability density according to 
( I83> at the end Ell (compare the similar situation for the 
variational ground-state wave function in Refs. I35ll36ll ). 



B. Anharmonic Oscillator (y > 0) 

The variational procedure described in the last section 
is now applied for determining the time evolution of the 
probability density in case of the nonlinear drift coefficient 
(0 with y > 0. By doing so, the optimization of the variational 
parameter K is performed for each value of the variable x at 
each time t. The result of such a calculation with the param- 
eters g — = 1,^ = = 0. 1 and /j — is shown in Figure 
|6l for the time interval < f < 20. Figure |6l (a) depicts the 
optimal variational parameter K which is the largest solution 
from ( 18 1> . For small times, the optimal variational parameter 
reveals a -dependence. For large values of f , the variational 
parameter becomes independent of x which corresponds to 
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the expectation that the probability density converges towards 
the stationary one. In Figure |S](b) the respective variational 
results for the distribution are compared with the cumulant 
expansion. The latter shows insofar a wrong behavior as it 
has two maxima whereas the numerical solution has only 
one. The variationally determined probability density is 



depicted by dots which correspond to the optimal variational 
parameters in Figure |5] (a). We observe on this scale nearly 
no deviation between the variational and the numerical 
distribution, thus already the first-order variational results are 
quite satisfactory. 




FIG. 6: Time evolution of probability density from variational optimization for g = D = 1,y=ct = 0.1, and jj = Q. (a) Largest variational pa- 
rameter K determined from <81> . colors code different times corresponding to the distributions shown in (b). (b) Dots correspond to variational 
parameters of (a) and coincide on this scale with numerical solutions of FP equation as represented by the lines through the dots. Cumulant 
expansions are shown as gray areas. At the front the stationary distribution iPstat(-^) and the corresponding potential 'I'(x) are depicted, (c) 
Distance <84> between variational and cumulant expansion from numerical solution of FP equation. 



In order to quantify the quality of our approximation, we 
introduce the distance between two distributions vix^t) and 
2'(x,f) at time t according to 

\,}if)^2j Fix,t)-'Pix,t)\dx. (84) 

If both distributions are normalized and positive, the maxi- 
mum value for the distance ^ (f ) is 1 and corresponds to the 
case that the distributions have no overlap. However, if they 
coincide for all x, the distance Aj, ^ (f ) vanishes. Thus small 
values of A^, -(f) indicate that both distributions are nearly 
identical. In Figure |6l (c) we compare the time evolution of 
the distance (jS^ between the variational distribution and the 
cumulant expansion from the numerical solution of the corre- 
sponding FP equation, respectively. We observe that the vari- 
ational optimization leads for all times to better results than 
the cumulant expansion, both being indistinguishable for very 
small and large times, as expected. 

C. Double Well (y<0) 

Now we discuss the more complicated problem of a non- 
linear drift coefficient (|3} with y < 0. The corresponding po- 



tential i2H has the form of a double well, i.e. it decreases 
harmonically for small x and becomes positive again for large 
X, so that a stationary solution exists (see the front of Figure 
Elb)). Strictly speaking, the cumulant expansion developed 
in Section ITV El makes no sense for y < as the unperturbed 
system g ^ does not have a normalizable solution. This 
problem is reflected in the time evolution of the cumulant ex- 
pansion coefficients P„.jt,/(f) which are listed in Table IIII and 
depicted in Figure |8] for the order n = I. In contrast to the 
case y > in ( I75> . the coefficient pi,4,o(f ) vanishes for y < 
so that the cumulant expansion (I69ll-( I7U does not even lead 
to the correct stationary solution j23L i2bi up to first order in 
g. In the limit t °° the only nonvanishing expansion coeffi- 
cients Pi,<:,/(f ) read 

y>0: limJ,o,(0^ 2y(D-ya') ' /l^J'^^^) = 2^^ 



11 




FIG. 7: Time evolution of probability density from variational optimization for g = 10, D = I, y = —1, a = 0.1 and /j = 0. (a) Largest 
variational parameter K determined from <81> . colors code different times corresponding to the distributions shown in (b). (b) Dots correspond 
to variational parameters of (a) and coincide on this scale with numerical solutions of FP equation as represented by the lines through the 
dots. Cumulant expansions are shown as gray areas. At the front the stationary distribution 2'stat(^) and the con'esponding potential 4>(x) are 
depicted, (c) Distance <84> between variational and cumulant expansion from numerical solution of FP equation. 




FIG. 8: Time evolution of the expansion coefficients ^i,kj{t) from 
Table IIforO= 1, Y= -1, A* = 0, and o = 0.5. 



The results of the first-order variational calculation of the 
probability density are summarized in Figure^Jfor the param- 
eters g = 10, D = 1, Y = — 1, o = 0.1, and /j ~ 0. Due to 
the strong nonlinearity, we could determine the optimal varia- 
tional parameter K from solving ( I81> for all x and t as shown 
in Figure0(a). As expected the cumulant expansion diverges 
for larger times t as illustrated in Figure (b). Despite of 
this the variational distribution lies precisely on top of the nu- 
merical solutions of the FP equation. This impressive result 
is also documented in Figure0c) where the cumulant expan- 
sion shows for increasing time t no overlap with the numerical 
solution, whereas the distance between the variational and the 
numerical distribution decreases. 

Even more difficult is the case of a weak nonlinearity where 
the two minima of the double well are more pronounced. 
Therefore, the variational calculation has also been performed 
for the parameter values g ~ 0.1, D — 1,Y=— l,o = 0.1, and 




FIG. 9: Zeros of 3pC (x,f;K)/3K with continuous solutions for t < 
fcrit = 1-618 and disconnected branches beyond. 



/J = 0. For small times one obtains again a continuous opti- 
mal variational parameter Kgpi(x,t) from solving ( I81> for all x. 
However, there exists a critical time fcrit ~ 1 -6 beyond which 
condition ( 18 U has different solution branches depending on x. 
In general such a critical point exists, if the equations: 




{-^— -^crit ^c-rit ^l^— l^ei-it } 







(86) 



have a solution {xcritifcritii^crit}- For the case at hand we find 
-fcrit = ±2.267, fcrit = 1-618 and Kcrit = —1.326- The corre- 
sponding surface of zeros and the critical points are depicted 
in Figure|9] As it remains unclear how these branches should 
be combined for evaluating the probability density, we resort 
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to zeros of higher derivatives (I82> that are continuous for all 
times. For small values of t we find one, two and three con- 
tinuous zeros for the first, second and third derivative, respec- 
tively, as shown in Figure^] for f = 1.2. The smallest zero 
of the second and third derivative reach a critical point at 
fcrit = 1.5284 and fcrit = 1.4588. The larger zeros have con- 
tinuous solutions for all values of t, as shown in FigurelTTIfor 
f = 1.7. 

Figure [21 shows the distance between variational and cu- 
mulant expansion from numerical solution of FP equation for 
these different branches of zeros. Apparently there is no easy 
choice of the right branch of zeros, which gives good results 
at all times. While Kopt,3b gives good results at small times, 
it fails to approach the stationary solution. The largest zero 
K«pt.3c of the third derivative on the other hand is unusable for 
small f , but gives good results for later times. 



Figure [O] shows the distance ( I87> for the branches of zeros 

















' fSffJI'""" 










-6 -4 


-2 2 


4 6 



FIG. 10: Zeros of <82t for ra = {1,2,3} at ? = 1.2, before reaching 
any critical time. 
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FIG. 12: Distance <84> between variational and cumulant expansion 
from numerical solution of FP equation for different branches of ze- 
ros. The smallest zero of each derivative reaches a critical time at 
/R^ 1.6. 
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FIG. 13: Distance <87> between moments of distributions deter- 
mined from variationally from different branches of zeros of the third 
derivative. 



FIG. 11: Zeros of {82} for m = {1,2,3} at f = 1.7. The smallest 
branch of zeros for each derivative has reached a critical point and 
dissected in disconnected branches. 

Assuming we had no knowledge of the numerical solu- 
tion, we need to find a way to switch the solution branch 
from Kopt,3b to Kopt,3c- In order to determine a suitable time 
to change the branch of zeros from Kopt(jc,f) to Kop((x,f), we 
consider the deviation of the moments of the corresponding 
distributions 



/_+r |x" [P(,^,f;Kopt) -f (x,f;<p,)] I dx 
J^::\x"P,mix,t)\dx 



(87) 



Kopt,3b(x,f) and Kopt,3c(^j0 for the first three even moments 
n —2,4,6. We find, that the distributions are in good agree- 
ment at f w 2.6, so we choose to combine the solution for 
Kopt.3b(-Y,r) for t < 2.6 with the solution for Kont3c(^,0 for 
t > 2.6. The combined result is shown in Figure^^ The dis- 
tance (I84> between variational expansion and numerical solu- 
tion of the FP equation for this case exhibits a small kink at 
t w 2.6 due to the change in the branch of zeros. Furthermore, 
in comparison with the other cases in Figures|6]andQthis dis- 
tance is relatively large which underlines that this is, indeed, 
a difficult variational problem. Note, however, that the com- 
bined solution succeeds in approaching the stationary solution 
for large times. 
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FIG. 14: Time evolution of probability density from variational optimization for g = 0.1, O = 1, Y = ~1. = 0.1, and /j = 0. (a) Variational 
parameter K determined from <82> for m = 3, colors code different times corresponding to the distributions shown in (b). For t < 2.6 we 
selected the branch of zeros Koptjb, whereas for t > 2.6 the branch Koptjc was used. Outlined dots are interpolated zeros (see text), (b) 
Dots correspond to variational parameters of (a) and coincide on this scale with numerical solutions of FP equation as represented by the lines 
through the dots. Cumulant expansions are shown as gray areas. At the front the stationary distribution 2'stat(^) and the corresponding potential 
<I>(x) are depicted, (c) Distance <84t between variational and cumulant expansion from numerical solution of FP equation. 



We remark that the variational approach of Ref. 12311 is related 
to ours. In contrast to our method, one obtains there in case 
of the difficult double well problem ^ = 0.1,K=— la unique 
solution of the extremal condition j8U for all x and t. How- 
ever, the resulting probability density shows for larger times t 
significant deviations from our, and from numerical solutions 
of the FP equation. 



D. Higher Orders 

High-order variational calculations have been performed 
for the double well with the parameters g~ 10, Z) — l,y=— 1 
in case of an initially Gaussian-distributed probability density 
peaked at the origin, i.e. o = 0.1,/j = for t = 0.23. This 
time was chosen due to its large distance between the varia- 
tional result and the numerical solution in order to reduce pos- 
sible errors in the numerical solution. The order of magnitude 
of the systematic error of the numerical solution can be esti- 
mated by comparing the numerical solution of the harmonic 
problem, e.g. g with the exact solution that is available 
for that case. We find that the error of the numerical solu- 
tion is about 10^^, which is smaller that the pointsize used in 
Figure [21 The first three variational orders, shown in Figure 
[Ts] converge exponentially to the numerical solution of the FP 
equation. 



ditive noise which is characterized by the nonlinear drift coef- 

10" 



10" 



10" 



1 



FIG. 15: Distance <84> between numerical solution of FP equation 
and the first three variational calculations forg = 10, Y=^l.£' = 1. 
o = 0.1,^ = Oandr = 0.23. 



ficient (|3jl. A comparison with numerical results shows an ex- 
ponential convergence of our variational resummation method 
with respect to the order We hope that VPT will turn out 
to be useful also for other applications in Markov theory as, 
for instance, the calculation of Kramer rates (see the re- 
cent variational calculation tunneling amplitudes from weak- 
coupling expansions in Ref. 0]), the treatment of stochastic 
resonance I3«l . or the investigation of Brownian motors |39]. 



VI. SUMMARY 
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